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Introduction 


In  this  papei  we  discuss  the  test  of  fit  that  a random  sample  comes 
from  the  extreme  value  distribution 

(l)  F(x)  = exp[-exp{-(x~£)/0}]  , < x < . 

This  is  the  most  important  of  the  three  distributions  which  arise  as  the 
asymptotic  distribution  of  the  largest  value,  suitably  normalized,  of  a 
sample  taken  from  any  of  a wide  class  of  distributions  (see  e.g.  Johnson 
and  Kotz  (1970),  p.  272),  It  arises  also  in  applied  research; 
the  present  work  was  motivated  by  a question  arising  in  biology.  Suppose 
the  sample  is  in  ascending  order  x^  <_  x^  <_  <_  * * • <_  x^  , and  the  null 

hypothesis  is  formally  stated  as  : the  sample  comes  from  distribution 

(l),  with  possibly  one  or  both  of  the  location  and  scale  parameters  £ , 

0 unknown.  Following  Stephens  (197^ , 1976)  we  distinguish  four  situa- 
tions : 

Case  0 : Both  £)  and  0 are  known,  so  that  F(x)  is  completely 

specified. 

Case  1;  0 known,  5 to  be  estimated. 

Case  2:  £ known,  0 to  be  estimated. 

Case  3:  £ , 0 both  unknown,  and  to  be  estimated. 

We  suppose  the  parameters  will  be  estimated  by  maximum  likelihood 
from  the  given  example;  the  estimates,  for  Case  3,  are  given  by  the  equa- 
tions (Johnson  and  Kotz  (1970 ),  p.  283): 

0 = exp(-x^/0)]/[£  exp(-x  /§)] 


(2) 


and 


(3)  t = -0  ln[Z^  exp(-x  /0)/n]  . 

Equation  (2)  is  solved  iteratively,  and  then  (3)  can  he  solved.  In  Case  1, 

A . A 

0 is  known;  then  £ is  given  by  (3)  with  0 replacing  6 . In  Case  2, 

£ is  known;  suppose  then  that  y = x_^  - E,  ; 0 is  given  by  solving 

(4)  0 = {Ey  - Ey  exp(-y, /0 )}/n  . 

<]  J J 3 J 

2.  The  Goodness-of-Fit  Tests 

The  tests  discussed  'below  are  based  on  EDF  statistics,  measuring  the 

discrepancy  between  the  empirical  distribution  function  and  the  theoretical 

distribution  (l);  in  Cases  1,  2,  and  3,  the  m.l.  estimates  are  inserted 

for  unknown  parameters  in  (l).  The  statistics  discussed  are  those  usually 

called  W2,  and  A2  . Asymptotic  theory  for  the  statistics  will  be 

given  in  Section  3;  this  is  based  on  work  of  Anderson  and  Darling  (1952), 

for  Case  0,  and  on  papers  by  Darling  ( 1955 ) » Sukhatme  (1972),  Durbin  (1973) 

and  Stephens  (1976)  for  situations  where  parameters  must  be  estimated. 

The  reader  is  referred  to  these  papers  for  the  theory  behind  the  methods 

which  follow;  in  particular,  we  follow  closely  the  general  lines  of 

Stephens  (1976).  In  that  paper,  for  example,  will  be  found  the  defini- 
? 2 2 

tions  of  W , U and  A ; here  we  shall  give  only  the  practical  steps 
in  making  a test  of  Hq  . These  are: 

(a)  Calculate  z ^ - F(x^)  where  F(x)  is  as  given  in  (l)  with  the 
appropriate  m.l.  estimates  inserted  for  unknown  parameters  in  Cases  1,  2, 
or  3.  Recall  that  the  x_^  are  supposed  in  ascending  order,  giving  z 
also  in  ascending  order.  Let  z be  the  mean  of  the  z ^ . 
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(b) 


Calculate  the  test  statistic  desired: 


o n 

w2  = l 


i=l 


(z. 

1 


2i-l'2 
2n  1 


+ 


1 

12n 


U2  = W2  - n(z  - 


A 


-[£. (2i-l){ln  z.  + ln(l“Z  . )}]/n 
l l n+l-i 


n . 


(c)  Refer  to  Table  1,  first  calculating  the  modified  statistic  and 
then  comparing  these  with  the  (upper  tail)  points  given  in  the  table, 
for  the  appropriate  case.  For  example,  if  a sample  of  size  10,  Case  3, 
gave  W2  - 0.119  » the  modified  W2  is,  say,  W*  =0.119(1-0. 2/J  ) = 0.1265; 

this  is  just  significant  at  the  level. 


3.  Asymptotic  Theory  of  the  Tests 

The  following  is  a summary,  taken  from  the  papers  referenced  above, 
of  the  steps  needed  to  calculate  the  asymptotic  distributions  of  W2 , U2 
or  A2  in  the  various  cases . 

For  Case  0,  when  F(x)  is  completely  specified,  it  is  well-known 

that  z = F(x)  gives  a random  variable  z which  is  uniformly  distributed 

between  0 and  1 ; further,  if  z^  = F(x  ) , i = l,2,...,n  , and  if 

F (z)  is  the  EDF  of  the  z.  sample,  then  y (z)  = /n{F  (z)  - z}  tends 
n inn 

asymptotically  to  a Gaussian  process  y(z),  0 < z < l , with 

E(y(z) ) = 0 for  all  z , and 

(5)  E(y(s)y(t))  = pQ(s,t)  = s - st  , 0 £ s 1 t £ 1 > 

where  E is  the  expectation  operator. 
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Unknown  parameters.  Suppose  now  that  the  (continuous)  distribution 
under  test  is  F(x)  (not  necessarily  (l)),  with  density  f(x)  , con- 
taining k parameters  0^ , 0O , . . . , 8^  ; these  are  unknown  and  will  be 
estimated  by  maximum  likelihood.  Suppose  Z is  the  matrix  with  entries 
Z. . given  by 


(6) 


= -E 


/ 3^  In  f(x)*\ 

\ 36i  36J  j 


9 


i s 0 ~ 1 »2, . . . ,k  5 


and  let  E be  its  inverse.  Further,  define  s = F(x)  , and 


(T) 


gjU) 


3s 

30.  * 


Let  u(s)  be  a k-vector  whose  i-th  component  is  g^(s)  . Then,  under 
appropriate  regularity  conditions  (see  e.g.  Durbin  Cl973)),  y^(z)  now 

tends  to  a Gaussian  process  y(z)  with  mean  zero  as  before,  and  with 
covariance 


(8) 


p(s,t)  = pQ(s,t)  - u’(s)Eu(t)  . 


The 
ditions , 
Anderson 
Gaussian 
function 
dropped ) 


statistic 


W2 


is 


0 


y^(z)  dz  and,  again  under  regularity  con- 


its asymptotic  distribution  will  be  that  of 


y (z)  dz  . 


and  Darling  (.1952)  and  Darling  (1955)  have  shown  how  another 
process  may  be  constructed,  with  mean  0 and  given  covariance 
p(s,t)  and  the  distribution  (the  word  asymptotic  will  be 
of  W2  is  calculated  from  the  new  process.  We  must  first  solve 


(9) 


f(x)  = X 


rl 

p(x,y)  f(y)  dy 

0 


k 


for  eigenfunctions  i\(y)  and  corresponding  eigenvectors  X.  . If  D(A) 
is  the  Fredholm  determinant  associated  with  (9),  the  characteristic  func- 
tion of  the  distribution  of  will  be  {D(2it )}-~/2  . distribution 

of  will  then  be  the  same  as  that  of 


(10) 


s = y v . /x. 

u 11 
i=l 


where  the  w,  are  independent  X-^  variables, 
distribution  are 


The  cumulants  of  the 


(11) 


K = 2d"1(j-l)l 


p.  (s,s)  ds 


where  p (s,t)  = p(s,t)  , and,  for  j >_  2 , p.(s,t)  is  defined  by 


Pj(s,t)  = 


rl 


p^_1(s,u)  p(u,t)  du 


The  cumulants  can  also  be  found  from  the  representation  (10) 


(12) 


K. 


= 2'j"1(j-l)l  l (1/A.  )3  . 

i=l  1 


2 2 

Statistics  U and  A . These  statistics  are  respectively  functionals 

of  y (z)_  y (z)  and  of  y (z)  w(z)  where  w(z)  = {z  (l-z )}-1/2  ; Y (z)-Y  (z) 
n - n n n n 

these  processes  tend  asymptotically  to  Gaussian  processes  with  covariances 

p(s,t)  , which,  for  a given  case,  can  be  found  from  the  corresponding 

covariance  for  ; then  Equation  (9)  must  again  be  solved,  and  the  sum 

2 2 

S in  (10),  with  the  new  , will  give  the  distribution  of  U or  A 
The  characteristic  function  comes  from  the  Fredholm  determinant  of  (9)  as 
before . 
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Thus  the  practical  problem  is  to  find  p(s,t)  for  different  cases, 
then  to  solve  (9)  for  X , and  then  to  use  (l2)  to  approximate  the 
distribution  (assuming,  as  is  nearly  always  so  except  for  Case  0,  that 
the  characteristic  function  cannot  be  inverted). 

4 . Asymptotic  Results  for  the  Extreme  Value  Distribution 

The  above  results  apply  to  any  suitably  regular  distribution.  For 
the  distribution  (l)  considered  here,  we  find,  after  some  algebra. 


2 

where  Y = .57712  is  Euler's  constant  and  c is 

r"(3)  + r"(i)  - 2r"(2)  - 1 = f2/6  + (y-i)2  = 1.82368  , 

giving  c = 1.350437  ; F"(x)  is  the  second  derivative  of  the  Gamma 

A A 

function  T(x)  . The  asymptotic  variances  of  £ (Case  l)  and  0 (Case  2), 

—1  2 —1 

when  0=1,  are  then  n and  (c  n)  , respectively.  The  functions 
g(s)  , with  subscript  to  indicate  the  case,  are 

gi(s)  = - e-  ’ 

_ / \ _ S In  s [-In (-In  s)3 

§2 ' S • Q 9 

and  these  give  covariance  functions  p (s,t)  for  Case  j,  j = 1,2  : 

d 

(i4)  Pj(s,t)  = pQ(s,t)  - 4>jCs)  4>j(t) 
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■with 


4>1  ( s ) = s In  s 

and 

<j>2(s)  = s In  s [-In (-In  s)]/l. 350^37  . 

2 2 -1 

For  Case  3 we  first  define  , cr^  , and  p by  -writing  Z = X 

as 


A A g 

The  large-sample  variances  of  £ and  0 , in  Case  3,  are  then  cr^/n  and 
2 

Og/n  respectively,  and  the  correlation  between  them  is  p . Inversion  of 
Z gives  (for  0=1),  = 1.10867  , = 0.60793  and  p = 0.313  . 

Hote  that  the  variances  are  different  from  those  in  Cases  1 and  2 because 
of  the  presence  of  correlation. 

Since  £ is  positive  definite,  it  can  be  written  £ = BB’  where  B 
is  upper  triangular;  an  example  given  by  Sukhatme  (1972)  is 


(16) 


If  we  define  v(s)  = B'u(s)  , with  two  components  i|^(s)  , ^g(s)  , the 

covariance  for  Case  3 becomes,  from  (8), 


(17)  P3(s,t)  = PQ(s,t)  - ^(s)  ^(t)  - l|>2(s)  Tp2(t ) . 
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Another  possible  matrix  decomposition  is  Z = CC'  where 


A A 

If  there  were  no  asymptotic  correlation  between  E,  and  0 in  Case 
3,  p0(s,t)  would  take  the  form  of  (17 ) , but  with  cj)^(»)  replacing 

V*)  • 


5 • Calculation  of  Weights 

With  covariance  p (s,t)  known  for  Case  j , the  next  problem  is  to 

J 

solve  (9)  for  the  weights  A^  . This  is  done  as  follows. 

Let  0 < \ < A2  < A3  •“  weiS^s  an^  f.(x)  the  associated 

normalized  eigenfunctions,  for  Case  0.  (For  W2  , these  are 
A.  = l/7r2i2  and  f . (x)  = /2  sin  rrix  .)  Set  D_(A)  = II  (A-A.)  . Expand 


4>l(x)  = £ a^^Cx)  , so  that  = (jj^^Cx)  fi(x)  dx  ; similarly  let 


b.  = 

l 


<j>2(x)  f (x)  dx  . Let 


Sa(X)  ' 1 + d-  • 

1 


b2 

sb(x)  = 1 + x J CTT  • 

1 


<wx)  = 
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In  general,  the  A_^  of  Cases  1 and  2 are  then  given  "by  setting  the 
Fredholm  determinant  to  zero;  this  implies,  for  Case  1,  solving 
Dq(X)  S&(X)  = 0 , and  for  Case  2,  Dq(X)  S^CX)  = 0 . For  Case  3, 
suppose  ip  (x)  , ip2(x)  replace  <J>  (x)  , 4>2(x)  in  the  definitions  of 

a . and  b . above , and  define 

T(A)  = SaU)  Sb(\)  - {Sab(X)}2  ; 

solutions  for  Case  3 are  then  given  by  Dq(X)  T(A)  = 0 . 

Weights  for  W2.  For  W2,  the  A^  of  Case  0,  which  we  call  the 

standard  weights,  each  occur  only  once.  In,  say.  Case  1,  if  a value  i 

exists  so  that  a_.  = 0 , then  A = A^  would  be  a solution  of 

D (A)  S (X)  = 0 ; otherwise  ( X— X . ) in  D-JA)  cancels  ( X— X . ) in  the 
□ a 10  i 

term  in  S (X)  and  X = X.  is  not  a solution.  For  distribution  (l) 
a l 

the  a^  are  never  zero  for  W2,  so  the  A^  are  then  given  only  by 
S (X)  = 0 . Similarly,  for  Case  2 we  solve  S,  (X)  = 0 , and  for  Case  3, 

8#  D 

T(A)  = 0 (see  Stephens  (l9?6)  for  examples  where  a^  = 0). 

2 2 2 

Weights  for  A and  U . For  A a similar  situation  exists;  the 
functions  f^(x)  for  Case  0 are  Pb(2x-l)  , where  P^(t)  are  Ferrer 
associated  Legendre  functions,  and  the  standard  A^  = i(i+l)  . No  coeffi- 
cient a.  or  b.  is  zero,  and  solutions  for  X.  for  Cases  1,  2,  and  3 
ii  i * 

are  given  by  ( X ) = 0 , S^CA)  = 0 > and  T(A)  = 0 , respectively. 

2 

For  U the  discussion  is  more  complicated.  The  roots  A^  of 

2 2 

D0(A)  = 0 are  double  roots,  given  by  X^  = Utt  i , and  the  corresponding 
eigenfunctions  are  f^(x)  = ^2  sin  2tt1x  and  f^(x)  = /2  cos  2irix.  Suppose  a 
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and  a*  , to . and  to*  are  the  coefficients  obtained  using  f.(x)  and 
11  i 1 

f*(x)  respectively.  Then  S (A)  becomes 
1 & 


a2  a*2 

8-(x>  “1  + xl?wr/r: + < • 

,x  11  1 


and 


K.  * 
a.  to . 
i i 


sab(x)  * M ? 13717  + 1 137X7 

,1  1 1 1 


S^(X)  is  defined  similarly  to  S^CA)  . In  Case  1,  although  no  coefficient 

a.  , a.  , to.  or  to.  is  zero,  one  factor  A - A.  in  D„(A)  cancels 
ill  i i 0 

denominators  in  S (A)  , tout  the  other  factor  remains,  so  A = A.  is  a 
a i 

solution  of  D (A)  S (A)  = 0 ; other  solutions  are  given  directly  toy 

U 3, 

S (A)  = 0 . Similarly,  for  Cases  2 and  3,  the  standard  weights  occur  once 
in  addition  to  those  given  respectively  toy  S^A)  = 0 and  T(.A)  = 0 . 

As  a check  on  calculations,  the  weights  for  Case  3,  for  all  three 
statistics,  were  found  using  the  c^(s)  and  cf^Cs)  given  toy  tooth  trans- 
formations B and  C of  Section  2;  these,  of  course,  give  different  coefficients 
* * 

a^ ,a^ ,tou  ,to.  tout  the  same  final  distribution  for  S in  (10). 


6.  Direct  Calculation  of  Means  and  Variances 

The  means  and  variances  of  the  various  distributions  can  toe  found 

directly  from  (ll).  For  example,  for  the  mean  y for  Case  j , (ll)  gives 

J 


(l8)  y^  = 


P1(s,s)ds  = 


Pq(s ,s )ds  - 


cj;1(s)ds  = yQ  - 


(j>2(s)ds 


The  mean  drops  from  its  Case  0 value  (y^  =1/6)  toy  the  last  integral 
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in  (18),  which  we  will  call  A^  ; this  can  sometimes  be  directly  cal- 
culated. For  example,  for  W2  Case  1,  using  integration  by  parts. 


Similarly, 


Let  s = e 


A1  ” 


<j>1Cs)  ds 


for  Case  2,  let  b be 


A_ 


cj)2(s)  ds  = b 


-z 


then 


fl 

0 

fl 

■*0 


s2ln2s  ds  = .07^07  . 
(l.350^37 )~2  ; then 
s2ln2s{ln(-ln  s)}2  ds  . 


A 


2 


= b 


1*00 

-3z  2.  2 . 

e z In  z dz  . 

•'0 


Integrals  of  this  type  can  easily  be  evaluated  by  use  of  appropriate 
substitutions  and  the  identities 


r(m)  = 
r’(m)  = 
r"(m)  = 


m-1 


,00 


m-1 

y 


Jo 


m-1 

y 


ay  » 

ID  y dy  , 
ln2y  dy  . 


For  Ag  above,  substitute  u = 3z  ; 


A -i 
2 27 


e~U  u2(ln  u - In  3)^  du 


0 


= b{r"(3)  - 2 In  3 r*(3)  + In  ‘3  r(3)}/27 


= 0.01750  . 
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For  W~  Case  3,  using  the  notation  of  Section  4,  after  Equation  (l6), 
we  have 


j-1 

• 0 
(•1 

•'0 

rl 


♦*(-) 


ds  + 


ds 


v'  (s ) v(s ) ds 
u'  (s ) £ u(s ) ds 


J0 

“ih  * 2pai02A12  + °2A2 


where  we  define 


rl 

4>  (s)  4)  ( s ) ds  . 
0 X 


Let 


c he 


1/1.350437  ; then 


s2lni's  {-ln(-ln  s ) }ds 

J0 


c 


rco 

e"3zz2  In  z dz  = 0.0096445  . 

•*0 


Thus  finally  = 0.1083  . The  three  means  for  Cases  1,  2,  3 are  then 
y^  = 0.0926  3 y^  = 0.1492  , y^  = 0.0584  . For  U2  the  calculations 
are  on  similar  lines  hut  are  somewhat  more  complicated;  the  results  are 
yQ  = 1/12  , \i±  = 0.0718  , y2  = 0.0683  , y3  = 0.0559  . For  A2  the 

integrals  are  intractable  and  have  been  calculated  numerically  to  giye 

yQ  = 1 , y1  = 0.5959  , y2  = 0.8619  , y3  = 0.3869  . 

In  principle,  variances  could  be  calculated  from  (ll),  as  was  done 
in  Stephens  (1976)  for  the  case  of  the  normal  distribution,  but  here  the 
integrals  were  too  complicated.  It  is  important,  however,  to  have  the 
exact  means,  as  (12)  usually  converges  too  slowly  to  give  them  correctly. 
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T. 


Calculation  of  Percentage  Points 


Wien  the  A^  are  found,  and  the  exact  means,  the  percentage  points 
of  S in  (10)  can  be  found  by  a modification  of  Imhof's  method,  given 
by  Durbin  and  Knott  (1972).  Alternatively,  the  first  four  cumulants  can 
be  found  from  (l2)  and  Pearson  curves  fitted  to  the  data.  Imhof's  method 
can  be  made  very  accurate  for  a finite  sum  in  (lO),  though  it  is  expen- 
sive in  computer  time;  in  adapting  it  for  an  infinite  sum  an  element  of 
approximation  is  introduced.  Thus  both  techniques  give  approximate 
percentage  points.  The  Imhof  method  gives  more  accurate  points  in  the 
lower  tail  (see  Solomon  and  Stephens  (1975)  for  a fuller  discussion),  but 
this  is  not,  of  course,  the  tail  which  would  generally  be  used  in  goodness 
of  fit  work.  Both  techniques  will  depend  on  the  accuracy  of  the  A^  , 
in  turn  dependent  on  the  accuracy  of  the  numerical  integrations  in  a^ 
and  b . . 

i 

Both  methods  were  used  here  to  find  the  percentage  points  given  in  Table  1; 
they  were  in  agreement  to  the  accuracy  given.  Points  for  Case  0 are  included  to 
show  how  much  the  points  drop  in  the  other  cases  (and  to  show,  therefore,  how 
important  it  is  to  use  the  correct  points). 

8.  Monte  Carlo  Results  for  Finite  n 

Monte  Carlo  studies  were  made  to  determine  the  percentage  points  of 
the  various  statistics,  for  sample  sizes  n = 10,  20,  and  50  ; 5,000 

samples  were  used  for  each  -case.  Previous  experience  had  suggested 
(Stephens  (l97^))  that  convergence  to  the  asymptotic  points  would  be 
rapid,  and  a plot  of  percentage  points  against  l/n  proved  this  to  be  so. 
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The  Monte  Carlo  points  were  used  to  calculate  the  modified  forms  given 
in  Table  1;  for  further  details  on  how  these  are  found,  see  Stephens 
(1970,  1976). 

Comment . It  is  interesting  that  the  asymptotic  process  in  Case  1 
above  has  exactly  the  same  covariance  (p  (s,t)  in  (l4))  as  in  a test  for 
the  exponential  distribution  F (x ) = 1 - exp(-x/0)  , with  0 unknown; 
see  Case  1+  of  Stephens  (197b).  The  asymptotic  distributions  of  all  the 
EDF  statistics  are  therefore  the  same  in  Case  1 above  as  in  the  exponen- 
tial test. 
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